#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
活化能计算器
基于Arrhenius方程计算活化能
参考: 1-s2.0-S0959652622013968-main.pdf中的计算方法
"""

try:
    import pandas as pd
    import numpy as np
    import openpyxl
    from openpyxl.styles import Font, PatternFill, Alignment
except ImportError:
    print("请安装必要的依赖包: pip install pandas numpy openpyxl")
    exit(1)

import warnings
warnings.filterwarnings('ignore')

class ActivationEnergyCalculator:
    def __init__(self, excel_path):
        """初始化活化能计算器"""
        self.excel_path = excel_path
        self.raw_data = None
        self.results = {}
        
    def load_data(self):
        """加载实验数据"""
        try:
            # 尝试读取第一个sheet
            self.raw_data = pd.read_excel(self.excel_path, sheet_name=0)
            print(f"成功加载数据，共{len(self.raw_data)}行")
            print(f"数据列: {list(self.raw_data.columns)}")
            return True
        except Exception as e:
            print(f"加载数据失败: {e}")
            return False
    
    def detect_columns(self):
        """自动检测温度、时间和转化率列"""
        columns = self.raw_data.columns
        
        # 寻找可能的列名
        temp_cols = [col for col in columns if any(word in str(col).lower() 
                    for word in ['temp', '温度', 'temperature', 't', '°c'])]
        time_cols = [col for col in columns if any(word in str(col).lower() 
                    for word in ['time', '时间', 'min', 't'])]
        conversion_cols = [col for col in columns if any(word in str(col).lower() 
                          for word in ['conversion', '转化率', 'α', 'alpha', 'conv'])]
        
        # 如果没有找到明确的列，使用前几列
        if not temp_cols:
            temp_cols = [columns[0]] if len(columns) > 0 else []
        if not time_cols:
            time_cols = [columns[1]] if len(columns) > 1 else []
        if not conversion_cols:
            conversion_cols = [columns[2]] if len(columns) > 2 else []
            
        return {
            'temperature': temp_cols[0] if temp_cols else columns[0],
            'time': time_cols[0] if time_cols else columns[1],
            'conversion': conversion_cols[0] if conversion_cols else columns[2]
        }
    
    def calculate_activation_energy_friedman(self, data, temp_col, conv_col):
        """使用Friedman方法计算活化能"""
        # 获取唯一温度
        unique_temps = data[temp_col].unique()
        
        results = []
        
        for temp in unique_temps:
            temp_data = data[data[temp_col] == temp].copy()
            temp_data = temp_data.sort_values(by=conv_col)
            
            # 计算转化率导数（简化方法）
            if len(temp_data) > 2:
                temp_data['dα/dt'] = temp_data[conv_col].diff() / temp_data.index.to_series().diff()
                temp_data['dα/dt'] = temp_data['dα/dt'].fillna(0)
            else:
                continue
            
            # 计算ln(dα/dt) vs 1/T
            temp_k = temp + 273.15  # 转换为开尔文
            temp_data = temp_data[temp_data['dα/dt'] > 0].copy()
            
            if len(temp_data) > 2:
                temp_data['ln(dα/dt)'] = np.log(temp_data['dα/dt'])
                temp_data['1/T'] = 1.0 / temp_k
                
                # 简单线性回归
                x = temp_data['1/T'].values
                y = temp_data['ln(dα/dt)'].values
                slope, intercept = self._linear_regression(x, y)
                activation_energy = -slope * 8.314  # R = 8.314 J/(mol·K)
                
                results.append({
                    'temperature': temp,
                    'activation_energy_kJ/mol': activation_energy / 1000,
                    'slope': slope,
                    'intercept': intercept,
                    'data_points': len(temp_data)
                })
        
        return pd.DataFrame(results)
    
    def _linear_regression(self, x, y):
        """简单线性回归"""
        n = len(x)
        sum_x = np.sum(x)
        sum_y = np.sum(y)
        sum_xy = np.sum(x * y)
        sum_x2 = np.sum(x * x)
        
        slope = (n * sum_xy - sum_x * sum_y) / (n * sum_x2 - sum_x * sum_x)
        intercept = (sum_y - slope * sum_x) / n
        
        return slope, intercept

    def calculate_activation_energy_kissinger(self, data, temp_col, conv_col, time_col):
        """使用Kissinger方法计算活化能"""
        # 计算每个温度下的最大反应速率对应的温度
        unique_temps = data[temp_col].unique()
        
        peak_data = []
        for temp in unique_temps:
            temp_data = data[data[temp_col] == temp].copy()
            temp_data = temp_data.sort_values(by=time_col)
            
            if len(temp_data) < 3:
                continue
                
            # 计算反应速率（简化方法）
            temp_data['dα/dt'] = temp_data[conv_col].diff() / temp_data[time_col].diff()
            temp_data['dα/dt'] = temp_data['dα/dt'].fillna(0)
            
            # 找到最大反应速率对应的温度
            max_idx = temp_data['dα/dt'].idxmax()
            if not pd.isna(max_idx):
                peak_temp = temp_data.loc[max_idx, temp_col]
                peak_data.append({
                    'temperature': temp,
                    'peak_temperature': peak_temp,
                    'peak_temperature_K': peak_temp + 273.15
                })
        
        if len(peak_data) < 2:
            return None
            
        peak_df = pd.DataFrame(peak_data)
        
        # 计算ln(β/Tp²) vs 1/Tp
        peak_df['ln(β/Tp²)'] = np.log(1.0 / (peak_df['peak_temperature_K'] ** 2))
        peak_df['1/Tp'] = 1.0 / peak_df['peak_temperature_K']
        
        # 线性回归
        x = peak_df['1/Tp'].values
        y = peak_df['ln(β/Tp²)'].values
        slope, intercept = self._linear_regression(x, y)
        activation_energy = -slope * 8.314 / 1000  # kJ/mol
        
        return {
            'activation_energy_kJ/mol': activation_energy,
            'slope': slope,
            'intercept': intercept,
            'peak_data': peak_df
        }
    
    def calculate_integral_method(self, data, temp_col, conv_col):
        """使用积分法计算活化能"""
        unique_temps = data[temp_col].unique()
        
        results = []
        
        for temp in unique_temps:
            temp_data = data[data[temp_col] == temp].copy()
            temp_data = temp_data.sort_values(by=conv_col)
            
            if len(temp_data) < 3:
                continue
                
            # 计算积分
            temp_k = temp + 273.15
            
            # 使用Coats-Redfern方法
            temp_data = temp_data[temp_data[conv_col] > 0.1].copy()
            if len(temp_data) > 2:
                temp_data['g(α)'] = -np.log(1 - temp_data[conv_col])
                temp_data['ln(g(α)/T²)'] = np.log(temp_data['g(α)'] / (temp_k ** 2))
                temp_data['1/T'] = 1.0 / temp_k
                
                # 线性回归
                x = temp_data['1/T'].values
                y = temp_data['ln(g(α)/T²)'].values
                slope, intercept = self._linear_regression(x, y)
                activation_energy = -slope * 8.314  # J/mol
                
                results.append({
                    'temperature': temp,
                    'activation_energy_kJ/mol': activation_energy / 1000,
                    'slope': slope,
                    'intercept': intercept,
                    'data_points': len(temp_data)
                })
        
        return pd.DataFrame(results)
    
    def process_all_methods(self):
        """处理所有计算方法"""
        if self.raw_data is None:
            print("请先加载数据")
            return False
        
        columns = self.detect_columns()
        print(f"检测到的列: {columns}")
        
        # 清理数据
        data = self.raw_data.copy()
        data = data.dropna(subset=[columns['temperature'], columns['conversion']])
        
        # 转换温度列为数值
        try:
            data[columns['temperature']] = pd.to_numeric(data[columns['temperature']], errors='coerce')
            data[columns['conversion']] = pd.to_numeric(data[columns['conversion']], errors='coerce')
            data = data.dropna()
        except:
            print("数据转换失败，使用原始数据")
        
        # 计算各种方法的活化能
        print("正在计算Friedman方法...")
        friedman_results = self.calculate_activation_energy_friedman(
            data, columns['temperature'], columns['conversion']
        )
        
        print("正在计算积分法...")
        integral_results = self.calculate_integral_method(
            data, columns['temperature'], columns['conversion']
        )
        
        print("正在计算Kissinger方法...")
        kissinger_results = self.calculate_activation_energy_kissinger(
            data, columns['temperature'], columns['conversion'], columns['time']
        )
        
        # 存储结果
        self.results = {
            'friedman': friedman_results,
            'integral': integral_results,
            'kissinger': kissinger_results,
            'columns': columns,
            'raw_data': data
        }
        
        return True
    
    def save_results(self):
        """保存计算结果到Excel"""
        try:
            # 加载现有工作簿或创建新的
            try:
                workbook = openpyxl.load_workbook(self.excel_path)
            except:
                workbook = openpyxl.Workbook()
                if 'Sheet' in workbook.sheetnames:
                    del workbook['Sheet']
            
            # 创建新的工作表
            if 'Activation_Results' in workbook.sheetnames:
                del workbook['Activation_Results']
            
            results_sheet = workbook.create_sheet('Activation_Results')
            
            # 写入标题
            results_sheet['A1'] = '活化能计算结果汇总'
            results_sheet['A1'].font = Font(bold=True, size=14)
            results_sheet.merge_cells('A1:E1')
            results_sheet['A1'].alignment = Alignment(horizontal='center')
            
            # 写入原始数据信息
            row = 3
            results_sheet[f'A{row}'] = '原始数据信息:'
            results_sheet[f'A{row}'].font = Font(bold=True)
            row += 1
            
            results_sheet[f'A{row}'] = f'数据行数: {len(self.results["raw_data"])}'
            results_sheet[f'A{row+1}'] = f'温度列: {self.results["columns"]["temperature"]}'
            results_sheet[f'A{row+2}'] = f'时间列: {self.results["columns"]["time"]}'
            results_sheet[f'A{row+3}'] = f'转化率列: {self.results["columns"]["conversion"]}'
            
            row += 5
            
            # Friedman方法结果
            if not self.results['friedman'].empty:
                results_sheet[f'A{row}'] = 'Friedman方法结果:'
                results_sheet[f'A{row}'].font = Font(bold=True)
                row += 1
                
                for col_idx, col_name in enumerate(self.results['friedman'].columns, 1):
                    results_sheet.cell(row=row, column=col_idx, value=col_name)
                    results_sheet.cell(row=row, column=col_idx).font = Font(bold=True)
                
                row += 1
                for _, data_row in self.results['friedman'].iterrows():
                    for col_idx, value in enumerate(data_row, 1):
                        results_sheet.cell(row=row, column=col_idx, value=value)
                    row += 1
                
                # 计算平均活化能
                avg_e_friedman = self.results['friedman']['activation_energy_kJ/mol'].mean()
                results_sheet[f'A{row}'] = f'Friedman方法平均活化能: {avg_e_friedman:.2f} kJ/mol'
                results_sheet[f'A{row}'].font = Font(bold=True)
                row += 2
            
            # 积分法结果
            if not self.results['integral'].empty:
                results_sheet[f'A{row}'] = '积分法结果:'
                results_sheet[f'A{row}'].font = Font(bold=True)
                row += 1
                
                for col_idx, col_name in enumerate(self.results['integral'].columns, 1):
                    results_sheet.cell(row=row, column=col_idx, value=col_name)
                    results_sheet.cell(row=row, column=col_idx).font = Font(bold=True)
                
                row += 1
                for _, data_row in self.results['integral'].iterrows():
                    for col_idx, value in enumerate(data_row, 1):
                        results_sheet.cell(row=row, column=col_idx, value=value)
                    row += 1
                
                avg_e_integral = self.results['integral']['activation_energy_kJ/mol'].mean()
                results_sheet[f'A{row}'] = f'积分法平均活化能: {avg_e_integral:.2f} kJ/mol'
                results_sheet[f'A{row}'].font = Font(bold=True)
                row += 2
            
            # Kissinger方法结果
            if self.results['kissinger']:
                results_sheet[f'A{row}'] = 'Kissinger方法结果:'
                results_sheet[f'A{row}'].font = Font(bold=True)
                row += 1
                
                results_sheet[f'A{row}'] = f'活化能: {self.results["kissinger"]["activation_energy_kJ/mol"]:.2f} kJ/mol'
                results_sheet[f'A{row}'].font = Font(bold=True)
                row += 1
                results_sheet[f'A{row}'] = f'斜率: {self.results["kissinger"]["slope"]:.4f}'
                results_sheet[f'A{row+1}'] = f'截距: {self.results["kissinger"]["intercept"]:.4f}'
                row += 3
            
            # 最终活化能汇总
            results_sheet[f'A{row}'] = '活化能汇总:'
            results_sheet[f'A{row}'].font = Font(bold=True, color='FF0000')
            row += 1
            
            if not self.results['friedman'].empty:
                avg_e_friedman = self.results['friedman']['activation_energy_kJ/mol'].mean()
                results_sheet[f'A{row}'] = f'Friedman方法: {avg_e_friedman:.2f} kJ/mol'
                row += 1
            
            if not self.results['integral'].empty:
                avg_e_integral = self.results['integral']['activation_energy_kJ/mol'].mean()
                results_sheet[f'A{row}'] = f'积分法: {avg_e_integral:.2f} kJ/mol'
                row += 1
            
            if self.results['kissinger']:
                results_sheet[f'A{row}'] = f'Kissinger方法: {self.results["kissinger"]["activation_energy_kJ/mol"]:.2f} kJ/mol'
                row += 1
            
            # 计算综合活化能
            all_values = []
            if not self.results['friedman'].empty:
                all_values.extend(self.results['friedman']['activation_energy_kJ/mol'].values)
            if not self.results['integral'].empty:
                all_values.extend(self.results['integral']['activation_energy_kJ/mol'].values)
            if self.results['kissinger']:
                all_values.append(self.results['kissinger']['activation_energy_kJ/mol'])
            
            if all_values:
                final_e = np.mean(all_values)
                results_sheet[f'A{row}'] = f'最终活化能: {final_e:.2f} kJ/mol'
                results_sheet[f'A{row}'].font = Font(bold=True, size=12, color='FF0000')
            
            # 调整列宽
            for column in results_sheet.columns:
                max_length = 0
                column_letter = column[0].column_letter
                for cell in column:
                    try:
                        if len(str(cell.value)) > max_length:
                            max_length = len(str(cell.value))
                    except:
                        pass
                adjusted_width = min(max_length + 2, 50)
                results_sheet.column_dimensions[column_letter].width = adjusted_width
            
            # 保存工作簿
            workbook.save(self.excel_path)
            print(f"结果已保存到 {self.excel_path}")
            
            return True
            
        except Exception as e:
            print(f"保存结果失败: {e}")
            return False
    
    def generate_report(self):
        """生成计算报告"""
        report = []
        report.append("=" * 50)
        report.append("活化能计算报告")
        report.append("=" * 50)
        
        if self.results:
            report.append(f"数据文件: {self.excel_path}")
            report.append(f"数据行数: {len(self.results['raw_data'])}")
            report.append("")
            
            if not self.results['friedman'].empty:
                avg_e = self.results['friedman']['activation_energy_kJ/mol'].mean()
                report.append(f"Friedman方法平均活化能: {avg_e:.2f} kJ/mol")
            
            if not self.results['integral'].empty:
                avg_e = self.results['integral']['activation_energy_kJ/mol'].mean()
                report.append(f"积分法平均活化能: {avg_e:.2f} kJ/mol")
            
            if self.results['kissinger']:
                report.append(f"Kissinger方法活化能: {self.results['kissinger']['activation_energy_kJ/mol']:.2f} kJ/mol")
            
            # 计算综合活化能
            all_values = []
            if not self.results['friedman'].empty:
                all_values.extend(self.results['friedman']['activation_energy_kJ/mol'].values)
            if not self.results['integral'].empty:
                all_values.extend(self.results['integral']['activation_energy_kJ/mol'].values)
            if self.results['kissinger']:
                all_values.append(self.results['kissinger']['activation_energy_kJ/mol'])
            
            if all_values:
                final_e = np.mean(all_values)
                report.append("")
                report.append(f"最终活化能: {final_e:.2f} kJ/mol")
        
        return "\n".join(report)

def main():
    """主函数"""
    excel_path = "/Users/rowanelsewhere/Desktop/crystal/ExpDat_1-7.xlsx"
    
    print("开始活化能计算...")
    calculator = ActivationEnergyCalculator(excel_path)
    
    if calculator.load_data():
        if calculator.process_all_methods():
            if calculator.save_results():
                print("计算完成！")
                print(calculator.generate_report())
            else:
                print("保存结果失败")
        else:
            print("计算过程失败")
    else:
        print("数据加载失败")

if __name__ == "__main__":
    main()